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The nonequilibrium dynamics of diffusion-mediated plasticity and creep in materials subjected to 
constant load at high homologous temperatures is studied atomistically using Phase Field Crystal 
(PFC) methods. Creep stress and grain size exponents obtained for nanopolycrystalline systems, 
m ~ 1.02 and p ~ 1.98, respectively, closely match those expected for idealized diffusional Nabarro- 
Herring creep. These exponents are observed in the presence of significant stress-assisted diffusive 
grain boundary migration, indicating that Nabarro-Herring creep and stress-assisted boundary mi¬ 
gration contribute in the same manner to the macroscopic constitutive relation. When plastic 
response is dislocation-mediated, power law stress exponents inferred from dislocation climb rates 
are found to increase monotonically from m ~ 3, as expected for generic climb-mediated natural 
creep, to m ~ 5.8 as the dislocation density pd is increased beyond typical experimental values. 
Stress exponents m > 3 directly measured from simulations that include dislocation nucleation, 
climb, glide, and annihilation are attributed primarily to these large pd effects. Extrapolation to 
lower pd suggests that m ~ 4 — 4.5 should be obtained from our PFC description at typical exper¬ 
imental pd values, which is consistent with expectations for power law creep via mixed climb and 
glide. The anomalously large stress exponents observed in our atomistic simulations at large pd 
may nonetheless be relevant to systems in which comparable densities are obtained locally within 
heterogeneous defect domains such as dislocation cell walls or tangles. 

PACS numbers: 61.72.Bb, 61.72.Lk, 62.20.F- 81.40.Lm 


The tendency of a solid material to gradually and ir¬ 
reversibly deform or even flow under low loads at high 
temperatures is termed creep deformatiorf 1 Low loads 
and high temperatures in this context are a < a y and 
T > T m /2, respectively, where a y is the yield stress and 
T m is the equilibrium melting temperature of the mate¬ 
rial. This type of slow plasticity not only alters material 
microstructure, shape, and properties over extended time 
periods, but is also a primary cause of mechanical failure 
in materials such as gas turbine blades that operate in 
high temperature load bearing environments 1 3 . 

The plastic flow that occurs during creep can involve 
conservative defect evolution mechanisms, such as dis¬ 
location glide and grain boundary sliding, but is gen¬ 
erally facilitated by thermally activated defect motion 
along local stress and/or chemical potential gradients, 
and is therefore inherently diffusive in nature. Vacancy 
diffusion in particular tends to be a central facilitator 
of deformation, either directly during diffusional creep 
or indirectly during dislocation or power law creep, as 
discussed further in the following. Moreover, it is the 
collective evolution of different defect populations over 
multiple length and time scales that leads to the observed 
diverse macroscopic phenomenology of creep. If the char¬ 
acteristic time scales associated with the diffusion of in¬ 
dividual defects are almost entirely inaccessible to most 
atomistic descriptions, then the length and time scales as¬ 
sociated with collective defect diffusion and macroscopic 
deformation are decidedly unreachable with conventional 
approaches. Details of the atomic-scale mechanisms by 


which diffusional (low a) and, in particular, power law 
(higher a) creep occur therefore remain matters of spec¬ 
ulation in many cases. There is clearly a need to better 
characterize these mechanisms, to understand how they 
collectively generate mesoscale dislocation pattering phe¬ 
nomena observed during creep, and to understand how 
these mesoscale structures connect to macroscopic creep 
phenomenology. 

The aims of the first part of the present study are 
to reproduce the macroscopic phenomenology of diffu¬ 
sional creep in nanopolycrystals from long time scale 
atomistic simulations, and to examine the consequences 
of stress-assisted grain boundary migration!? 11 on diffu¬ 
sional creep. The primary aim of the second part of this 
study is to begin addressing the question of whether such 
simulations can also provide direct connections between 
atomistic defect evolution mechanisms and dislocation- 
mediated power law creep phenomenology. Toward these 
ends we employ the phase field crystal (PFC) modeling 
approach. We begin with a brief overview of creep phe¬ 
nomenology and a discussion of simulation literature in 
Section |!| In Section [H) the PFC approach is introduced 
and a new method for conducting constant stress PFC 
creep simulations is outlined. Grain boundary mediated 
deformation in dislocation-free nanopolycrystals is then 
examined in Section [TIT| while dislocation-mediated creep 
is examined in Section |IV| Connections to power law 
creep phenomenology are discussed in Section m and 
conclusions are summarized in Section El 
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I. CREEP PHENOMENOLOGY 

The macroscopic plastic response of a material under¬ 
going creep deformation can generally be characterized 
by an equation of the form 1 3 



where e is total plastic strain, t is time, p is shear mod¬ 
ulus, d is grain size, and m and p are stress and grain 
size exponents, respectively. In laboratory experiments, 
a constant stress is typically applied and the resulting 
macroscopic deformation is measured as a function of 
time t. The total strain can be decomposed into elastic 
and plastic contributions, e(t) = e el + e pl (£). A rapid elas¬ 
tic deformation is thus observed, followed by an initial 
regime of non-steady plastic flow called primary creep. 
Plastic flow eventually reaches a steady-state during the 
so-called secondary creep regime, such that the strain 
rate is constant in time and Eq. 0 is applicable. 

Diffusional creep, characterized by a stress exponent 
m ~ 1, is often observed at low a in such experiments. 
At sufficiently high T, bulk vacancy diffusion is the dom¬ 
inant flow mechanism and a grain size exponent p ~ 2 
results (Nabarro-Herring creepP). At lower T, vacancy 
diffusion along grain boundaries is the dominant flow 
mechanism and a grain size exponent p ~ 3 results (Coble 
creepl 6 ). Previous atomistic studies have addressed as¬ 
pects of diffusional creep. Coble creep has been re¬ 
produced in molecular dynamics (MD) simulations of 
be d 12 ! 13 and fee 14 nanopoly crystals, for example. This 
is feasible since vacancy diffusion rates can be relatively 
rapid within loosely packed high angle grain boundaries. 
Signs of non-negligible Nabarro-Herring creep contribu¬ 
tions have been reported as the grain size is increasecP^, 
but a clear transition to dominant Nabarro-Herring creep 
has not been observed in MD simulations to our knowl¬ 
edge. 

Deviations from the idealized descriptions of Nabarro- 
Herring and Coble creep are known to occur in systems 
with non-trivial microstructures. For example, when all 
grains in a polycrystal have the same shape and size, 
the idealized models of diffusional creep predict a sym¬ 
metric and affine elongation of all grains. However, the 
presence of a distribution of grain sizes has been shown 
to lead inevitably to stress-assisted grain growth, the 
shrinkage of some grains and growth of others, during 
diffusional creep 7 . Thus, when grain boundaries are suf¬ 
ficiently mobile, idealized diffusional creep behavior can 
be significantly modified by the simultaneous effects of 
stress-assisted grain boundary migration. Consequences 
include non-trivial, non-affine evolution of grain mor¬ 
phology and the suppression of grain elongation for suffi¬ 
ciently large grain boundary mobilities and low stressed 3 . 
These effects have been examined with both meP and 
mesoscale simulation method^. 

At applied stresses higher than those at which diffu¬ 
sional creep is typically observed, significant dislocation 
motion begins to occur as vacancy fluxes around disloca¬ 
tion cores and rates of stress-assisted thermally activated 
motion increase. A transition from diffusional to m > 1 


or power law creep then occurs, such that m ~ 3 — 10, 
with the precise value depending on a variety of material 
and microstructural properties, many of which are not 
well understood. Simultaneous interactions involving va¬ 
cancies, dislocations, grain boundaries, and other defects 
can lead to complex behaviors such as climb-assisted ob¬ 
stacle bypasd^ on atomistic scales and collective dislo¬ 
cation patterning®] on mesoscopic scales. All manner of 
such processes may contribute to macroscopic power law 
creep rates. 

The simplest description of dislocation-mediated 
power law creep, m = 3 natural creepP, is arrived at 
by assuming a linear relation between the steady-state 
dislocation velocity and stress, v ss ~ cr, and by also 
employing a phenomenological relation between dislo¬ 
cation density pd and stress (Taylor’s relation), pd ~ 
a 2 . Orowan’s equation for plastic strain then specifies 
6 pl (£) = y ^2^ =1 {bi£iS x i) : where V is the system vol¬ 
ume, N is the number of mobile dislocations, bi is the 
magnitude of the Burgers vector of dislocation i, £i is 
its length, and Sxi is the distance that it has traversed 
at time t. The corresponding incremental plastic strain 
rate is e pl (£) = (M^)> where vi is the velocity 

of dislocation i at time t. Substituting Vi -» v ss ~ a 
and N/V = pd ~ & 2 gives e pl ~ a 3 or the natural creep 
exponent mn — 3. Larger stress exponents m 4 — 6 are 
often observed even in pure metals, and are generally be¬ 
lieved to be associated with significantly more complex 
processes than those underlying natural creep. 

Atomistically-informed kinetic Monte Carlo (kMC) 
simulations have been used to study dislocation climb 
and power law creep in highly deformed bee Fe 17 . Steady 
state climb velocities computed from such simulations 
were found to be highly nonlinear in stress, increasing 
as v ss ~ a q with q ~ 3 — 3.5 for large dislocation den¬ 
sities, pd ~ 10 15 /m 2 -10 17 /m 2 . This strong nonlinearity 
was found to arise from enhanced vacancy diffusion rates 
at large pd as well as from the assumed vacancy supersat¬ 
uration, applicable to highly deformed metals. By then 
applying Orowan’s equation, e pl = pdbv ss , and assum¬ 
ing that Taylor’s relation holds in such systems, pd ~ cr 2 , 
the authors estimated power law creep exponents directly 
from the measured climb velocities as m ~ 5 — 5.5. This 
result essentially describes a modified form of natural 
creep, wherein large pd and a steady supersaturation of 
vacancies increase the ‘inherent’ dislocation climb stress 
exponent from q = 1 to q > 3. 

Clouet 18 has connected the atomistic modeling of Ref. 
m to classical mesoscale descriptions of climb, and 
found that extrapolation of their results to more phys¬ 
ically relevant dislocation densities leads to stress expo¬ 
nents m ~ 3, as generally expected for climb-assisted 
natural creep in metals. The additional effects of glide, 
at minimum, must also apparently be considered if exper¬ 
imental stress exponents m > 4.5 are to be reproduced. 

Other mesoscale descriptions have been employed to 
study creep deformation. Recent discrete dislocation dy¬ 
namics (DDD) simulations that incorporate both glide 
and climb have reproduced responses consistent with 
both diffusional and power law creep in fee metalJ 19 ! It 
therefore appears that the length scales, dislocation den- 




3 


sities, and underlying mechanisms relevant to some forms 
of power law creep may be accessible with this type of 
coarse grained description. A difficulty lies in the ab¬ 
sence of atomistic detail, which necessitates the use of 
numerous explicit rules for dislocation interaction scenar¬ 
ios. Guidance for these rules often comes from atomistic 
simulations, but these are not always universally applica¬ 
ble or even well-defined, and only a subset of all possible 
interactions can realistically be considered. Thus, atom¬ 
istic simulations that permit basic characterization of 
diffusion-mediated defect evolution processes should con¬ 
tribute to better informed mesoscale simulations. Long 
time scale atomistic methodologies that can access meso¬ 
scopic length scales through large scale simulations or 
systematic coarse graining procedures may also eventu¬ 
ally provide fully self-contained multi-scale descriptions 
of creep and plasticity phenomena. 


II. PFC APPROACH 


PFC models contain atomistic detail but describe time 
scales in crystals near and beyond that of the character¬ 
istic vacancy diffusion time 2 ( 3“ Both conservative and 
nonconservative dislocation evolution processes are cap¬ 
tured for arbitrary crystal structure, orientation, mor¬ 
phology, and applied stres^^ 31 l PFC simulations of 
diffusion-accommodated plastic flow under creep condi¬ 
tions may therefore reveal previously inaccessible infor¬ 
mation about the atomistic mechanisms that control dif- 
fusional and power law creep. 

A general PFC free energy functional can be written 


F = 



1 2 W 3 / —| u A 

-n (r) - -n (r) + — n (r) 


1 

2 



dr dr 2 n(r)C 2 {\r - r 2 |)n(r 2 ). 


( 2 ) 


where F = F / (ksT pi), ks is Boltzmann’s constant, pi 
is a constant reference density, n(r) = p(r)/pi — 1 is the 
scaled time-averaged atomic density field, p(f) is the un¬ 
sealed time-averaged atomic number density field, w and 
u are free coefficients, and — f^l) is the two-point di¬ 

rect correlation function of the fluid, assumed isotropic. 
n(r) may assume nonzero average values uq. We uti¬ 
lize the structural or XPFC class of functionals 32 33 , for 
which the Fourier transformed correlation function can 
be written 


C 2 (k)i = e -(k-k i ) 2 /(2a 2 i ) e -a 2 T k 2 i /(2 Pi p i )^ ( 3 ) 


where i denotes a family of lattice planes at wavenum¬ 
ber ki. The constants a^, pi, and fa are the Gaussian 
width (which sets the elastic constants), planar atomic 
density, and number of planes, respectively, associated 
with the ith family of lattice planes, gt is a Debye- 
Waller-type temperature parameter that modulates the 
scattering intensity (S(k) = (1 — C 2 (&)) _1 where S(k) is 
the structure factor) due to the effect of atomic thermal 
vibrations. The envelope of all selected Gaussians i com¬ 
poses the final (%(&). A single reflection^ at k\ = 2 


is used here to produc e eq uilibrium bee structures with 
lattice constant a ~ y/2/3. 

The standard stochastic nonlinear diffusion equation is 
employed for n(r) dynamics, 


dn(f) 

dt 



+ V(r,t), 


(4) 


where t = Tt/pj is rescaled time (denoted as t in subse¬ 
quent sections), T is a mobility constant, and r](r,t) is a 
Gaussian stochastic noise variable with (r](f,t)) = 0 and 
( 7 ?(fi,ti)r?(f 2 ,t 2 )) = —2k B TV 2 S(fi -f 2 )S{ii -f 2 ). PFC 
studies of mechanical properties often employ an equa¬ 
tion of motion with an additional inertial or wave-like 
term, d 2 n/dP , to facilitate rapid elastic relaxations 26 . 
Since the diffusive component of the mechanical response 
is our primary interest in this work, and since we focus 
on the low stress and low steady-state strain rate regime 
in which inertial effects should be secondary, we employ 
the simpler diffusive form of Eq. ©• No qualitative dif¬ 
ferences in terms of the creep-type response outlined in 
the following sections were observed in comparative sim¬ 
ulations employing inertial dynamics. 

We further define M = 2k bT and use this variable 
in place of ksT to specify a scaled temperature. This 
description thus contains two explicit temperature pa¬ 
rameters, gt and M, as well as any implicit temperature 
dependences of the other parameters in Eqs. © and <©• 
For simplicity and following convention, it is assumed 
that any implicit T dependences in these latter scaled 
parameters can be neglected to lowest order. M and gt 
can thus be varied simultaneously to control the physical 
temperature T in the vicinity of the melting temperature 
(T m ), though either parameter has the same basic effect 
independently for our purposes. For simplicity, we there¬ 
fore control T by varying only M while holding gt fixed, 
unless noted otherwise. All simulations in this study were 
performed in 3D using a pseudo-spectral algorithm with 
semi-implicit time stepping and periodic boundary con¬ 
ditions. 

Within standard PFC descriptions, the explicit evo¬ 
lution of vacancies is interpreted as being coarse grained 
into the structure of n(f) and its diffusive evolution. One 
question that we wish to address is whether this coarse 
grained description of vacancy diffusion can reproduce 
Nabarro-Herring and Coble creep mechanisms, which in¬ 
volve only the strain-induced flow of vacancies. 

It is also currently unclear, in terms of power law creep, 
whether the length scales, dislocation densities, and col¬ 
lective defect kinetics relevant to m > 3 kinetics can be 
reached by PFC-type models. The case of natural creep, 
as outlined in Section [IJ provides an illustration of this is¬ 
sue. Since the kinetics of climb in PFC models have been 
shown to follow v ss ~ g for small dislocation densit}®, 
it stands to reason that a realization of Taylor’s relation 
would lead directly to m = 3 natural creep, or poten¬ 
tially m = 4.5 creep when glide is also considered 35 . It 
has not yet been demonstrated that the system sizes and 
dislocation creation-annihilation kinetics needed to self- 
consistently reproduce phenomena such as Taylor’s rela¬ 
tion can be simulated by PFC models. Though we do 
not attempt to simulate large 3D systems with multiple 
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slip systems in which Taylor’s relation generally emerges, 
analagous behaviors in quasi-2D systems are examined 
further in the following. 

We numerically investigate the response of bee crys¬ 
tals and polycrystals to creep-type deformation as fol¬ 
lows. As in typical creep experiments, deformation is 
applied under constant stress conditions, using a newly 
developed stress-controlled PFC simulation method. At 
each evolution iteration At, the grid spacing is homo¬ 
geneously varied in one or more directions i, depending 
on the particular choice of boundary constraints, until a 
specified stress a a is achieved. We consider quasi- 2 D sys¬ 
tem geometries with columnar grains aligned along the 
x axis and strain e applied along the y and z axes (see 
Fig. 0 , and further employ a fixed y-z area constraint 36 
with a xx = 0. A 2 and A y are simultaneously varied as 
A 2 = A 6 (l + e) and A y = A 6 /(l + e), where A 6 is the ini¬ 
tial grid spacing along i, until the condition df/de = a a 
is achieved, where / = F/V is the mean total free en¬ 
ergy. A x is then similarly varied until df/dA x = 0, 
though generally A x A 6 at all times for the columnar 
geometry employed. This method is therefore effectively 
equivalent to a constant volume pure shear deformation 
with fixed a a = cr yy + g zz . 

We have confirmed that the total strain generated ver¬ 
sus g zz for the case g xx = &yy = 0 reproduces the correct 
linear stress-strain relation of the perfect one-mode bee 
crystal for e < 0 . 02 . As expected, modest deviations 
from linearity are observed for e > 0.02 when nonlinear 
elastic effects become appreciable. We have also verified 
that the one-mode bee Poisson ratio v — 1/3 is realized 
at low g zz in these simulations, again with nonlinear ef¬ 
fects causing a modest gradual decrease to v ~ 0.31 by 
e = 0.08. 


III. DIFFUSIONAL CREEP AND 
STRESS-ASSISTED GRAIN BOUNDARY 
MIGRATION IN NANOPOLYCRYSTALS 

An idealized columnar grain structure with initially 
uniform grain shapes and sizes was used to examine dif- 
fusional creep in bee nanopolycrystals. Periodic systems 
with either 4 or 6 hexagonal grains per repeat unit were 
chosen, with various symmetric and asymmetric combi¬ 
nations of grain rotations and grain boundary tilt angles. 
High angle boundaries were examined in asymmetric 4 
grain systems with grain rotations of (0°, 45°, ±22.5°) or 
(±12.5°, ±37.5°). The former configuration contains four 
45° and eight 22.5° boundaries, while the latter contains 
two 15°, two 25°, and eight 40° boundaries, as shown in 
Figs. 0a) and [ 2 |c) , respectively. Symmetric 6 grain sys¬ 
tems with grain rotations of (0°, ±30°) or (45°, ±15°) 
were also examined. These configurations contain eigh¬ 
teen 30° boundaries, as shown in Figs. [ 2 ja) and 0 b). 
Low angle grain boundaries were examined using a simi¬ 
lar asymmetric 6 grain configuration with grain rotations 
of (0°, ±5°). This system contains twelve 5° boundaries 
and six 10° boundaries (see Fig.[5|a) ahead). The colum¬ 
nar axis is along x = [ 100 ] in all cases, and a system 
thickness L x = la was used in this direction (no changes 




FIG. 1: (Color online) Diffusional creep in bee nanopolycrys¬ 
tals with high angle grain boundaries, (a) Strain normalized 
system configurations for a a = 0.01 and d = 50.8a are shown 
at (Ot, e = 0.0), (62.5t, e = 0.031), (125t, e = 0.060), and 
(187.5t, e = 0.090). The system thickness is L x = la. Strain 
normalization is given by y —)> y( 1 + e) and z z/( 1 + e). 
For analysis and visualization purposes, local peaks in n(r), 
which represent the most probable atomic positions, are taken 
to correspond to atomic sites. Sites with bee coordination (at 
t = 0) are shown in pale green, those with irregular coordina¬ 
tion (grain boundary atoms) are shown in red {t = 0), green 
(t — 62.5), pale blue (t — 125), and gold (t — 187.5). (b) A 
representative set of creep curves for the configuration of (a). 
e pl vs. time is shown at various a a (color legend). The dashed 
black lines are linear fits. 


were observed for L x = 10 a). Grain sizes d = 25.4a, 
50.8a, 76.3a, 89.0a, 101.7a, 114.4a, 127.1a, and 152.5a 
were studied, where d is the diameter of the smallest 
circle that encloses the hexagon. Unless specified other¬ 
wise, model parameters w = 1.4, u = 1, no = 0, aq = 1, 
gt = 0.1, pi = 1, /?i = 8 , and M — 0 are used with 
A 6 ~ a /12 and At = 0 . 01 , where a = 0.81675. 

After equilibration of each system at g xx = & yy = 
&zz — 0 , simulations were conducted at a range of ga 
values, and the subsequent deformation e(t) = [L z (t) — 
L z (0)]/L z (0) was monitored as a function of time. Anal- 





















5 


t = 0 t = 300 

t = 600 

t = 900 

tb 

x^ z 

— °ZZ — *- 

SS | f 

CO 

O 

o 

0° 

G 

o 

CO 

1 

o 

O 

CO 

1 

30° 

0° 


t = 0 t = 500 t = 1000 t = 1500 



FIG. 2: (Color online) Morphology of diffusional creep 

in other bee nanopoly crystal configurations with high angle 
grain boundaries. Strain normalized system configurations 
are shown for: (a) The symmetric 6 grain cell with rotations 
of (0°, ±30°) and cr A = 0.01, d m 76.3a. Red: (0 1, e = 0.0), 
green: (300t, e = 0.044), pale blue: (600t, e = 0.090), gold: 
(900t, e = 0.135). (b) The symmetric 6 grain cell with ro¬ 
tations of (45°, ±15°) and ga — 0.01, d — 76.3 a. Red: 
(01 , e — 0.0), green: (500t, e = 0.108), pale blue: (lOOOt, 
e m 0.234), gold: (1500t, e — 0.384). (c) The asymmetric 4 
grain cell with rotations of (±12.5°, ±37.5°) and ga = 0.01, 
d = 101.7 a. Red: (Ot, e = 0.0), green: (450t, e = 0.055), pale 
blue: (900t, e = 0.111), gold: (1350*, e = 0.166). 


ysis of creep dynamics is restricted to stresses above cr_, 
the threshold for observable plasticity (<r_ ~ 0.0006) and 
below the threshold for dislocation nucleation from 
grain boundaries (cr + ~ 0.06). The absence of lattice 
dislocations or other defects in the grain interiors allows 
us to isolate and characterize grain boundary-mediated 
plasticity mechanisms without complications from col¬ 
lective dislocation processes, for example. Plastic flow 
in these grain geometries generally requires a significant 
amount of nonconservative grain boundary motion, i.e., 
vacancy-mediated migration and climb. 

Even though all grains are initially identically hexag¬ 
onal, their different orientations and boundary struc¬ 
tures create small asymmetries upon equilibration, which 
grow when stress is applied. The result is an initial 
delta function grain size distribution that broadens with 
time due to stress-assisted grain boundary migration 
or grain growth. The symmetric 6 grain geometry de¬ 
scribed above has been used in MD to avoid precisely 
this effect!^, but it will be shown in the following that the 
mobility of grain boundaries in our PFC simulations is 
large enough to significantly impact the measured steady- 
state creep rates. In discussing evolution of the grain size 
distribution, we distinguish between net grain growth, an 
increase in the average grain size, and differential growth, 
a broadening of the grain size distribution without any 
increase in average size. Net growth is observed only at 
very late times when entire grains are eliminated, after 
the steady-state creep regime has ended, and therefore 
does not influence our results. Differential growth does 
occur during the steady-state creep regime. 


A. High angle grain boundaries 

Representative creep curves obtained for the asymmet¬ 
ric 4 grain system with (0°, 45°, ±22.5°) rotations are 
shown in Fig. [ljb) for d = 50.8a. After an initial elastic 
strain e el and a brief regime of plastic flow onset (primary 
creep), an extended regime of linear steady-state creep 
(secondary creep) is observed for all cr_ < a a < cr + . The 
slope of the best linear fit to this regime gives the steady- 
state creep rate e ss for a given a a- The steady-state creep 
regime terminates in one of two ways. For a a < 0.06, dif¬ 
ferential grain growth proceeds until two grains eventu¬ 
ally consume the others at late times (net grain growth), 
leaving two 45° boundaries along y that no longer mi¬ 
grate in response to a a- Creep therefore ceases and all 
strain energy is absorbed elastically. For a a > 0.06, dis¬ 
locations eventually nucleate from the grain boundaries 
or homogeneously from the grain interiors at late times. 
This leads to unbounded plastic flow, a divergence in the 
strain rate, and loss of mechanical integrity. 

The extracted steady-state creep rates e ss for four val¬ 
ues of d are plotted as a function of a a in Fig. [3^a). The 
best fit slopes of these data for cr_ < a a < 0 + (dashed 
lines) indicate that the bare stress exponents are m = 
1.12 =b 0.02. Taking the effect of cr_ into account by sub¬ 
stituting ga — ct- for ga produces m = 1.02 ±0.02 (solid 
lines and inset). The dependence of e ss on d for 5 rep¬ 
resentative values of ga is plotted in Fig. |3|b). All data 
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FIG. 3: (Color online) Creep exponents of the configuration 
shown in Fig. [l] (a) Compilation of steady-state creep rates 
vs. a a at various grain sizes d. The dashed lines are fits to Eq. 
0 , the solid lines are fits to Eq. 0 with cr a replaced by a a — 
<7 -. The resulting creep rate stress exponents m indicated in 
the figure are those of the solid lines, while the corresponding 
values for the dashed lines are m — 1.11, 1.10, 1.14, and 
1.13, respectively. The inset shows the same data plotted vs. 
cJA — cr- with the same solid line fits, (b) Dependence of e ss on 
grain size d at various a a- Solid lines are fits to Eq. 0 with 
variable grain size exponent p. See Supplemental Material 37 
for animations of the 4 simulations at a a = 0.01. 
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FIG. 4: (Color online) (a) Effect of stochastic temperature 
M ~ ksT on grain boundary structure. Colors are as in Fig. 
[l|a). (b) Variation of the grain size exponent p with M. 


are examined further at the end of this section. We also 
note that grain rotations are observed in coordination 
with differential grain growth, and are a consequence of 
nonconservative boundary migration. 

Interestingly, we observe a transition from m ~ 1 and 
p ~ 2 behavior to m ~ 1 and p ~ 3 as the stochas¬ 
tic temperature M is raised to very near melting. The 
effects of M on grain boundary structure and the cor¬ 
responding grain size exponents p are shown in Fig. [4] 
There is a clear correlation between the increase in p 
from ^2 to ^3 and the effective liquification of the grain 
boundaries. The observed exponents are consistent with 
those of Coble creep, which is typically dominant only 
at relatively low T. Our findings therefore indicate that 
grain boundary premelting in the presence of low applied 
stresses can result in a re-entrant Coble creep regime at 
T sufficiently near the equilibrium melting point. 


sets are well fit by a grain size exponent p = 1.98 =b 0.05. 

Very similar values for the exponents m and p are ob¬ 
tained for the other high angle grain boundary configu¬ 
rations shown in Fig. [2] The morphology of the grain 
structures, on the other hand, does vary with the de¬ 
tails of the initial configuration. For example, the initial 
systems shown in Figs. [2|a) and 0b) are identical ex¬ 
cept for a rotation of all grains by 45° about the x axis, 
but their structures evolve quite differently. Whereas the 
two 0° grains in Fig. |2^a) grow nearly radially while the 
±30° grains shrink, all six grains in Fig. [2|b) elongate 
and coherently translate nearly uniformly along y. In the 
asymmetric 4 grain systems shown in Figs.[lja) an d§c), 
a relatively fixed 30° orientation is maintained along the 
boundaries that do not lie parallel to y. These bound¬ 
aries translate nearly uniformly along y until two of the 
four grains are eliminated. The causes of these differences 


B. Low angle grain boundaries 

The 6 grain low angle grain boundary configuration 
evolves as shown in Fig. [5ja). Since the Burgers vectors 
of the grain boundary dislocations can be clearly identi¬ 
fied, the direction of the driving force for each boundary 
under a zz and a yy is readily known (white arrows in Fig. 
|5ja)). Indeed, the motion of each boundary is consis¬ 
tent with that expected from its Burgers vector. Devia¬ 
tions from the exact direction expected are due to forces 
exerted by neighboring boundaries, as discussed further 
in the following subsection. Motion of the dislocations 
within boundaries aligned with y occurs purely by climb, 
while motion of all other boundary dislocations occurs 
via a mixture of climb and glide. The stress exponent 
m obtained during the regime of morphologies roughly 
corresponding to those shown in Fig. |5ja) is m 1.19 
(Fig. [5jb)). Eventually, boundaries begin to meet and 
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FIG. 5: (Color online) Diffusional creep in bee nanopolycrys¬ 
tals with only low angle grain boundaries, (a) Strain normal¬ 
ized system configurations for a a = 0.03 and d = 101.7a are 
shown at (Ot, e = 0.0), (40t, e = 0.0211), (lOOt, e = 0.0342), 
and (180t, e — 0.0499). The system thickness is L x — la. 
Sites with bee coordination (at t = 0) are shown in pale green, 
those with irregular coordination (dislocation core atoms) are 
shown in red (t = 0), green (t = 40), pale blue (t =s 100), and 
gold (t = 180). Burgers vectors of dislocations along each 
boundary are indicated with white _L symbols, along with the 
subsequent climb directions expected upon deformation, (b) 
Steady state creep rates vs. a a at d = 101.7a. The dashed 
line is a fit to Eq. 0 with m = 1.19. (c) Left: Creep curves 
for the data shown in (b). e pl vs. time is shown at various a a 
(color legend). The dashed black lines are linear fits. Right: 
Dependence of e ss on grain size d at a a = 0.01 and 0.03. Solid 
lines are fits to Eq. |l]) with grain size exponent p — 1.68. 


both the dislocation density and strain rate decrease at 
later times as dislocations annihilate. Figure |5|c) shows 
the grain size dependence of e ss , from which an exponent 
p ~ 1.68 is obtained. 


C. Discussion 

The creep exponents m ~ 1 and p ~ 2 obtained in 
the preceeding simulations are consistent with those of 
Nabarro-Herring creep, grain redistribution mediated by 


bulk vacancy diffusion. As discussed in Section [TJ ide¬ 
alized diffusional creep would be expected to generate a 
nearly affine deformation morphology, such that no sig¬ 
nificant structural changes would be observed upon nor¬ 
malization of the system configuration by strain (y —)> 
2/(1 + e), z z/(l + e)), no matter the initial grain ori¬ 
entations and patterns. This is clearly not the case in 
any of the present simulations. The grain deformation 
seen in Fig. [2jb) is the closest to affine, but some elon¬ 
gation along y is observed as well a significant amount 
of collective translation in this direction. During the pri¬ 
mary and steady-state deformation regimes, the summed 
grain boundary length and the number of grains within 
each system remain fixed, but the grain size distribu¬ 
tions broaden. Differential grain growth or stress-assisted 
grain boundary migration therefore clearly influences the 
measured creep rates. In this subsection, we attempt 
to characterize the relative degrees to which Nabarro- 
Herring creep and stress-assisted grain boundary migra¬ 
tion mediate plastic flow in these simulations, and to ex¬ 
plore connections, similarities, and differences between 
these two related mechanisms. 

The low angle grain boundary simulations are per¬ 
haps most easily understood from the viewpoint of stress- 
assisted migration. The expected climb direction of each 
dislocation/boundary, based on its Burgers vector and 
the direction of applied stress, is indicated by a white 
arrow in Fig. [5ja). The general agreement with the ob¬ 
served migration directions indicates that the applied 
stress directly induces dislocation climb and therefore 
boundary migration. The entire structure then evolves 
collectively in response to these climb forces, while the 
topological constraints of the boundary network transmit 
additional network forces to each boundary. For exam¬ 
ple, the dislocations within the upper left boundary of 
the central 0° grain migrate in the -z direction due to 
the climb force generated by a yy , but they also migrate 
in the +y direction due to the motion of the neighbor¬ 
ing vertical boundary in this direction. The net result 
is a mixed climb-glide migration nearly along the vector 
y — z. The two 10° vertical boundaries translate only 
along y by pure climb, as the forces due to neighboring 
boundaries balance along z. The four 5° vertical bound¬ 
aries climb similarly along y , but also experience a net 
torque about x due to the opposing migration directions 
of their neighbors along y and - y , respectively. 

By tracking the evolution of individual dislocations, 
we have quantified the amount of plastic strain relief due 
to dislocation motion. We find that this mode of strain 
relief accounts for nearly all of the applied system-wide 
strain, confirming that dislocation/boundary migration 
is the primary mechanism of plasticity in this case. 

Though high angle grain boundaries cannot be ana¬ 
lyzed in terms of individual dislocations, the relatively 
similar differential grain growth morphologies and creep 
exponents observed across low and high angle systems 
suggest that no qualitative change in plasticity mecha¬ 
nisms occurs with angle. The probability of conserva¬ 
tive grain boundary sliding should nonetheless become 
more significant as the boundary angles increase. Vi¬ 
sual inspection of n(r) during high angle boundary mi- 













gration reveals that density peaks are created within 
boundaries under tension (positive mass flux) and de¬ 
stroyed within boundaries under compression (negative 
mass flux), while no signs of significant boundary slid¬ 
ing have been observed. This flow of mass is consistent 
with the basic Nabarro-Herring process. It is therefore 
perhaps not surprising that the Nabarro-Herring expo¬ 
nents, m = 1 and p = 2, are observed even though the 
grain morphology does not remain affine during concur¬ 
rent stress-assisted migration. 

These findings suggest that plasticity in the present 
PFC simulations is mediated by the inherently coupled 
processes of vacancy diffusion and grain boundary mi¬ 
gration, in the following general manner. Applied stress 
generically induces vacancy flow out of grain boundaries 
under tension and into those under compression. These 
flows should in turn be capable of driving both the grain 
redistribution of diffusional creep and the nonconserva¬ 
tive defect evolution of stress-assisted boundary migra¬ 
tion. In our simulations, diffusion apparently occurs pri¬ 
marily through the grain interiors, in accord with ideal¬ 
ized Nabarro-Herring creep. Rather than simply deposit¬ 
ing onto non-tensile boundaries such that affine grain 
redistribution and elongation occurs, the vacancies di¬ 
rectly facilitate climb-mediated grain boundary migra¬ 
tion. This migration is possible because of the large 
nonconservative mobilities of PFC dislocations and grain 
boundaries relative to the inherent diffusive time scales of 
vacancies in this description. The subsequent directions 
of grain boundary migration are determined at the atom¬ 
istic level by the dislocation orientations and/or grain 
boundary topologies. These should generally exhibit no 
predisposition for affine morphologies, as is clear for the 
case of Fig. [5ja). The symmetry of the grain shapes 
is thus broken, a distribution of sizes emerges, and dif¬ 
ferential growth proceeds, driven by the stress-induced 
vacancy fluxes. The correspondence in terms of underly¬ 
ing mechanism suggests that this type of creep plasticity 
in the large boundary mobility regime may reasonably 
be expected to produce the same stress and grain size 
exponents as those of idealized Nabarro-Herring creep, 
though with very different effects on grain morphology, 
particularly for nanoscale grains. 

The non-affine morphologies appear to emerge here 
largely from atomistic effects, many aspects of which 
would therefore be difficult to predict without atomistic 
insight. First, variation of grain boundary structure with 
angle leads to breaking of the perfect hexagonal symme¬ 
try under zero stress. Second, the preferred and realized 
kinetics of boundary migration under stress are deter¬ 
mined by a number of potential effects with atomic-scale 
origins. For example, the stress dependence of the grain 
boundary mobilities and their preferred migration direc¬ 
tions are linked to grain boundary/dislocation structure. 
The driving force for migration may also be influenced 
by grain orientation if elastic anisotropy is sufficiently 
large. Further investigation is needed to determine which 
of these or other factors most influence the nature of this 
type of stress-assisted boundary migration, particularly 
for high angle boundaries, but over the range of parame¬ 
ters and deformation conditions investigated, our results 


indicate that it can be a significant contributor to creep 
flow in the regime of large boundary mobilities. 


IV. DISLOCATION OR POWER LAW CREEP 

To examine the transition to dislocation-mediated 
plastic flow and potentially power law creep, we must, at 
minimum, introduce dislocation sources into our simula¬ 
tions. Naturally emerging steady-state pd and e values in¬ 
volving primarily dislocation climb and glide would then 
signal a physically meaningful power law creep regime. In 
the following subsection, we begin by examining a simple 
controlled geometry in which the number of mobile dis¬ 
locations N does not depend significantly on cr^, and a 
stress exponent m = 1 is therefore expected (this is still 
technically dislocation, not diffusional, creep). It will be 
shown that m 1 can nonetheless arise due to addi¬ 
tional factors associated primarily with large dislocation 
densities. Other scenarios are then considered in which 
N does increase with a a and larger stress exponents are 
thus expected, as in the case of natural creep. Connec¬ 
tions are drawn between the subsequent power law stress 
exponents and the dislocation density effects identified in 
the constant N case. 


A. Dislocation nucleation, climb, and pile-up 

To begin examining the effects of dislocation climb on 
PFC creep rates, we first characterize a simple system 
geometry consisting of a perfect bee crystal with one cen¬ 
tral dislocation source, as shown in Fig. [6] The simula¬ 
tion cell and crystallographic directions are the same as 
those of the previous section, as is the system thickness 
of L x = la. The dislocation source is a cylindrical region 
of radius R = 10a in which n(f) is slaved to a homoge¬ 
neous penalty function? 7 . Though the model is capable of 
spontaneously generating dislocations from defects such 
as sessile dislocation loops and grain boundaries 3 ^, cylin¬ 
drical sources are input here “by hand” primarily for con¬ 
venience. They facilitate relatively controlled dislocation 
nucleation behavior at stresses below the homogeneous 
nucleation threshold and from sources small enough to 
introduce several within a single simulation cell. The 
cylinder acts as a rigid, uniform body similar to a large 
glass bead in a colloidal crystal or an incoherent inclusion 
in an atomic crystal. When a sufficiently large stress a a 
is applied as previously discussed, the cylinder becomes 
a site for heterogeneous nucleation of b = a(100) disloca¬ 
tion pairs or sets of orthogonal pairs (see Supplemental 
Material 37 ^ for animation). Since the stresses are axial, 
these dislocations nucleate and translate by climb alone. 
We can therefore characterize the creep response due to 
purely diffusion-mediated dislocation motion. In the case 
shown in Fig. [6j hard wall boundaries are also employed 
by slaving n(r) near the system edges to the perfect bee 
solution 27 . 

If the strain decomposition e(t) = e el + e pl (t) is em¬ 
ployed once more for the constant a a condition, then 
e pl (£) can be quantified using the Orowan equation when 
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FIG. 6: (Color online) Dislocation nucleation, climb, and pile-up in a simple dislocation creep scenario at <ja — 0.09 and 
L y m L z = 264 a. (a) 50 1, e = 0.0358, (b) 130t, e =s 0.0397, (c) 360t, e ==t 0.0544, (d) 480t, e = 0.0751. Interior atoms with 
bee coordination are shown in pale green, exterior hard wall atoms (also bee coordination) in pale blue, those with irregular 
coordination (dislocation core atoms) in red, and those with fee coordination in bright green. Red atoms are displayed with 
larger radii to highlight the dislocation positions. 


the only plastic strain relief mechanism is dislocation mo¬ 
tion, as noted previously. If N and ii are both constant 
and the average steady-state dislocation velocity can be 
written v ss ~ a q A , then e pl ~ For pure climb, q = 1 
is generally assumed and has been verified to hold in 
PFC simulation^!. Thus we expect m = q = 1 in con¬ 
stant stress experiments for which the dislocation density 
pd = N/V is small and constant and the overall stresses 
and dislocation velocities are low and linear, respectively. 
This should approximately be the case between nucle¬ 
ation of the first and second waves of dislocation pairs in 
large L y ,L z systems. Analyses and comparisons of ob¬ 
served nucleation rates and pile-up spacings with those 
predicted by continuum descriptions are outlined in Ap¬ 
pendices [A] and [Bj respectively. The agreement is good 
in both cases. 

Figure [7]^a) shows the measured dependence of m = q 
on system size or dislocation density pd. The exponent 
appears roughly to approach 1 as pd —» 0, though sys¬ 
tem size limitations prevent us from accessing very low 
pd values. For easily accessible p^, these simple dislo¬ 
cation creep simulations exhibit m > 2. The deviation 
from m = 1 appears to have two primary causes. The 


first, as already noted, is large p^. For a lattice con¬ 
stant of 0.33nm (as assumed in the following), the den¬ 
sities simulated correspond to pd — 10 14 /m 2 — 10 16 /m 2 , 
which is near or above the highest average values typ¬ 
ically observed in experiments (^10 12 /m 2 ). Our simu¬ 
lations therefore indicate that collective interactions be¬ 
tween defects at large pd can significantly accelerate plas¬ 
tic flow relative to that expected for isolated dislocations. 
A total pd < 10 12 /m 2 should therefore ideally be main¬ 
tained in such simulations, though large local densities 
pd ~ 10 12 /m 2 to 10 16 /m 2 may have direct physical rele¬ 
vance in some cases. 

The other source of m > 1 is associated with nonlin¬ 
ear acceleration of v ss at large as shown in Fig.^b). 
This acceleration correlates with the onset of significant 
nonlinear elasticity and can be avoided for a a < 0.05, 
but the existence of a finite stress barrier ajy for disloca¬ 
tion nucleation (ay - 0.035 in this case) also limits the 
explorable range at low a a- 

If we follow the analysis of Ref. m by applying 
Taylor’s relation directly to our measured climb veloc¬ 
ities, then the predicted creep rates would follow e ^ 

PdbVss ~ (j 2 +[ 1 - 8 ’ 3 - 8 ] ~ a [3.8,5.8] f or ^he va l ues of p d ex _ 
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FIG. 7: (Color online) (a) Stress exponent m for simple 
dislocation creep vs. dislocation density pd- Inset: Steady 
state strain rate data from which m values were determined 
(L = L y = L z ). (b) Steady-state dislocation velocity v ss vs 
applied stress a a for the first 6 dislocation pairs in the lowest 
Pd system examined. The red lines are fits to the first pair of 
nucleated dislocations that climb along d -z (see Fig. ©»)• 


amined. These values are comparable to those obtained 
from the atomistically-informed kMC simulations of Ref. 
H2, m ^ 5 — 5.5, which were conducted at dislocation 
densities near the high end of our pd range. The under¬ 
lying causes of these large exponents were identified in 
their study as nonlinear vacancy-dislocation core inter¬ 
actions and a non-negligible dependence of vacancy con¬ 
centration on stress. Such mechanisms are not explicitly 
considered in PFC descriptions, but nonlinearities near 
dislocation cores are naturally captured by n(r), the lo¬ 
cal peaks of which decrease in amplitude and increase 
in width as the core is approached. These variations 
signify local changes in vacancy concentration and mi¬ 
gration barriers, which appear to produce qualitatively 
similar nonlinearities in climb rates with increasing pd to 
those observed in the atomistic simulations of Ref. im. 

For sufficiently low pd (< 10 13 /ra 2 ), our PFC descrip¬ 
tion appears to converge toward the standard description 
of purely climb-mediated creep in which q = 1 and thus 
m = 3. The noted nonlinearities caused by interacting 
dislocation strain fields apparently become negligible at 
and below typical physically relevant dislocation densi¬ 
ties. This finding is in accord with Clouet’s analysid^ 
of the results of Ref. m wherein m ~ 3 is obtained 
for experimentally relevant conditions. It therefore ap¬ 
pears that much larger, less controlled system configu¬ 
rations in which dislocations may more realistically nu¬ 
cleate, climb, glide, annihilate, etc. must be examined 
to obtain a more meaningful description of power law 


creep. The required length scales are not yet accessible 
with PFC simulations, but we initiate steps in this direc¬ 
tion in the following section. We note nonetheless that 
the overall consistency of our results with those of atom¬ 
istic and mesoscale descriptions indicates that both the 
atomistic mechanisms and long-range interactions that 
underlie dislocation climb are qualitatively captured by 
PFC models. 


B. Power law creep 


We next examine somewhat less idealized microstruc¬ 
tures to determine whether more varied modes of dis¬ 
location creation, interaction, and annihilation can self- 
consistently generate power law creep behavior within 
atomistic PFC descriptions. As a first case, dislocation 
sources are symmetrically intr oduced into the hexagonal 
grain structures of Section III as shown in Fig.[8ja). This 
configuration produces a crossover from m ~ 1 diffusional 
creep for a a < <?n to m ~ 2.8 for a a > as shown 
in Fig.[8|b) (see Supplemental Material 37 for animation). 
Below ct/v, dislocations do not nucleate and the response 
is unchanged from that outlined in Section |III| for dif¬ 
fusional creep and stress-assisted grain boundary migra¬ 
tion. Above (T/v, dislocations nucleate from the sources 
and the response therefore transitions toward a disloca¬ 
tion mediated power-law-type regime. This change in 
strain relief mechanisms above cr at is clearly visible, with 
a crossover to larger but still nearly linear steady-state 
strain rates after dislocations become active. The mo¬ 
bile dislocation density becomes large during this regime, 
on the order of pd — 10 15 /m 2 , and increases roughly as 
Pd ~ The observed exponent m ~ 2.8 therefore 

seems to result primarily from the combined effects of 
this pd ~ <Ja dependence and the large pd and nonlinear 
v ss effects, as illustrated in Fig.[7| Glide is active in these 
systems, but does not appear to contribute significantly 
to the plastic flow in this idealized system geometry. 

A second system configuration containing 16 randomly 
positioned sources with L y = L z = 352a and 4 hard 
wall boundaries is shown in Fig. [9] (see Supplemental 
Material 37 ^ for animation). Such systems nucleate waves 
of dislocations when a a > ctn, and between early waves 
generally reach temporary steady-states in terms of p s J 
and e ss . Between the first and second waves of disloca¬ 
tion nucleation, for example, when pd — 3x 10 _15 /m 2 (at 
all stresses) we obtain m cx 4.2 (Fig. |9jb)). This value is 
larger than that expected from large pd and nonlinear v ss 
effects, which as quantified in Fig. [7] would be expected 
to produce m 3.0 =b 0.4. Similar simulations with 2 
rather than 4 hard wall boundaries produce m ~ 4.0. 
The excess m in these simulations may be due to the 
contribution of dislocation glide. 

If we extrapolate these results to lower pd where q = 1 
is expected to be realized (decrement measured m by 
~2), and further assume that Taylor’s relation holds (in¬ 
crement measured m by 2), a power law creep stress 
exponent of m ~ 4.0 — 4.2 is implied. This value ap¬ 
proaches but is smaller than typical experimental results 
for pure metals, m > 4.5, and is approximately consis- 
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FIG. 8: (Color online) Combined dislocation and 22.5°/45° 
grain boundary creep at d = 152.5 a with 1 central dislocation 
source per grain, (a) System representation at a a = 0.04 and 
t — 160, with colors as in Fig. [6] (b) Steady state strain rate 
data and stress exponents m. Inset: The corresponding creep 
curves at various a a (color legend). 


tent with previous mesoscopic descriptions that consider 
both climb and glide (m = 4.5) 3 ^. These findings, though 
not obtained directly at the system sizes and dislocation 
densities ulitimately sought, do indicate that physically 
meaningful behavior is being captured. 

As noted previously, the power law exponents obtained 
here at large pd may also have some limited direct phys¬ 
ical relevance based on the known spatial heterogeneity 
of pd in many systems. Even though average values of 
Pd > 10 -5 /a 2 ~ 10 14 /m 2 simulated here are not typ¬ 
ically encountered in real materials, significant spatial 
heterogeneity in Pd(f) is often observed. Dislocation clus¬ 
tering at cell walls, tangles, obstacles, and even low-to- 
medium angle grain boundaries, for example, have been 
estimated^ to produce local values of pd ~ 10 16 /m 2 , 
as large as any considered here. Though these hetero¬ 
geneous domains of very large pd are separated by do¬ 
mains of much lower pd (^ 10 13 /m 2 in Ref. [38]), they 
may nonetheless locally induce strong deviations from 
the q = 1 climb kinetics commonly assumed for indi¬ 
vidual dislocations, based on the simulations of Section 
IV A| Apparent power law creep behavior may therefore 
be observed in such systems, even in the absence of any 
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FIG. 9: (Color online) Dislocation creep at L y = L z — 352 a 
with 16 randomly positioned sources and 4 hard wall bound¬ 
aries. (a) System representation at a a = 0.08 and t — 65, 
with colors as in Fig. [6] (b) e ss data from which m was deter¬ 
mined. Inset: Creep curves at various a a for the configuration 
shown in (a). 

particularly complex defect reactions or mesoscale pat¬ 
terning. 


V. CONCLUSIONS 

We have studied the nonequilibrium kinetics of creep 
deformation and diffusional defect migration in atomistic 
crystalline and nanopolycrystalline systems using phase 
field crystal modeling. A method for conducting con¬ 
stant stress PFC simulations, the first such method to our 
knowledge, was developed and utilized throughout this 
study. The characteristic stress and grain size exponents 
quantified for symmetric nanopolycrystals, m ~ 1.02 
and p ~ 1.98, respectively, closely match those expected 
for idealized diffusional Nabarro-Herring creep, m = 1 
and pm 2. We find that a significant portion of the 
plastic flow in these systems is associated with non- 
affine grain boundary motion, suggesting that concurrent 
stress-assisted diffusive boundary migration does not nee- 
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essarily alter diffusional creep exponents in the regime 
of large boundary mobilities. Exponents consistent with 
diffusional Coble creep (m = 1, p = 3) were observed 
only in the presence of stochastic thermal noise ampli¬ 
tudes M large enough to induce grain boundary melting 
or premelting. This suggests that a re-entrant transition 
to Coble creep may occur in weakly stressed systems that 
exhibit significant premelting near T m . 


In simulations of dislocation-mediated plastic flow 
(Section |I V A ), nucleation rates and pile-up spacings 
were shown to be well described by predictions of con¬ 
tinuum elastic dislocation theory, Eqs. (A3) and (Bl), 
respectively. At numerically accessible dislocation den¬ 
sities, strong nonlinear interactions were found to pro¬ 
duce ‘inherent’ dislocation climb exponents 1.8 < q < 3.8 


and extrapolated creep stress exponents 3.8 < m < 5.8 
for pd ~ 10 14 /a 2 to 10 16 /a 2 . These results are in gen¬ 
eral agreement with previous atomistic 17 and mesoscale 
studies^, and further highlight the need to access lower 
dislocation densities in atomistic simulations. Extrapo¬ 
lated climb rates appear to approach q ~ 1 near typical 
experimental dislocation densities (pd ~ 10 12 /a 2 ), indi¬ 
cating that PFC models will reproduce climb-mediated 
natural creep behavior with m = 3 when Taylor’s relation 
holds. 


In systems with slightly less i dealize d grain and dislo¬ 
cation source stuctures (Section |IVB ), stress exponents 
m ~ 2.8 — 4.2 were directly measured, but these values 
cannot be straightforwardly compared to experimental 
results due to the known influence of large pd effects. 
Anomalously large stress exponents associated with the 
kinetic climb rate nonlinearities observed at large pd may 
nonetheless have relevance to power law creep in sys¬ 
tems that exhibit comparable dislocation densities locally 
within dislocation tangles and cell walls, for example. 
Extrapolation of the simulation results of Section |IVB| 
to physically relevant dislocation densities leads to pre¬ 
dicted exponents ofmc^ 4.0 —4.2, nearly consistent with 
mesoscale predictions for systems in which both climb 
and glide are active (m = 4.5), but still somewhat smaller 
than typical experimentally measured values in pure met¬ 
als (m > 4.5). 


The results of this work indicate that PFC-type ap¬ 
proaches are capable of capturing the essential physics 
of nonconservative defect evolution and dislocation- 
mediated creep at atomistic length scales and diffusive 
time scales. The system sizes required to describe typical 
experimentally-relevant dislocation densities, though not 
yet accessible, may be realizable within the next decade 
or potentially sooner with further development of compu¬ 
tationally efficient, coarse grained complex amplitude ex¬ 
pansions of PFC models^®*!. Potential points of shorter- 
term study include investigation of fully 3D systems with 
more realistic microstructures and a distribution of dislo¬ 
cation source activation stresses, as well as development 
of constant stress deformation methods that impose fewer 
constraints on the shape of the polycrystalline domain. 
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Appendix A: Dislocation Nucleation Rates 

As dislocation pairs sequentially nucleate in the sim¬ 
ulations of Section |IV A[ the dependence of pd on time 
and a a can be understood as follows. Upon nucleation 
of the first dislocation pair at £i, the local strain around 
the source is reduced by e zz = b/L z . For a a — cr at, this 
amount of strain must be reapplied to nucleate another 
pair, and the time required to do so, £2 — £ 1 , is equal to 
the time it takes the climbing dislocation pair to plasti¬ 
cally relieve this same amount of strain. For a constant 
steady-state climb velocity, 


t 2 -ti = 


L v 


LA l 2t, ss - 


(Al) 


At £2 there are 4 mobile dislocations, so the time required 
to nucleate the third pair is halved, and so on. In general, 


U ti — 1 — 


L v 


L z ifz 2(i - IK 


(A2) 


and the time required to generate i pairs of dislocations 
is 
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2c<JA z ' n 

n= 1 


(A3) 


where c ~ 6.59 is the constant prefactor given in Fig. 
0 b )- Equation (|A3|) is a harmonic series, which grows 


faster than logarithmically, implying that pd diverges 
faster than exponentially when annihilation and/or im¬ 
mobilization do/does not occm®l This would be the 
case, for example, in an infinitely large system. Since our 
simulations are finite and also employ hard wall bound¬ 
aries, dislocations become effectively immobile when they 
reach a wall. Nucleation of a new pair tends to occur 
roughly when the previous pair reaches the hard wall 
boundary, such that the number of mobile dislocations 
N remains nearly con stant with time. For this type of 
fixed N scenario, Eq. (A3) predicts that periodic waves 


of dislocation nucleation will occur every L y /(2v ss ) such 
that the total number of mobile plus immobile disloca¬ 
tions increases linearly with time at a rate proportional 
to (TA/Ly. These relations are indeed found to describe 
the first few waves of nucleation quite well for all system 
sizes studied with hard wall boundaries. 
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Appendix B: Dislocation Pile-Up Spacings 


Figure [To] shows the linear dislocation density p p (z) = 
fz-dz N( z ')dz'/ (2 dz) as a function of distance from the 
source for pile-ups such as that seen in Fig. (6^c). p p (z) 
is defined such that the number of dislocations between 
z — dz and z + dz is 2 dzp p {z). The solid line is a best 
fit to the prediction of continuum elasticity theory 44 ! for 
a double ID pile-up along z, 


Pp(z) 


2(1 — u)(ja z 

Jb y/{L z /2)* - 


(Bl) 


where z is the distance from the center of the pile-up. 
This equation with v = 1/2 (as imposed by the fixed 
area, constant a a algorithm), a a = 0.045, p, = 1, and 
b = 2/3 accurately describes the simulation results for 

z sufficiently far from the source. This further confirms 
that the long range elastic interactions between PFC dis¬ 
locations agree with continuum elastic descriptions, such 
that non-trivial multi-dislocation, multi-obstacle interac¬ 
tions are captured as well. 



FIG. 10: (Color online) Comparison of dislocation pile-up 
spacing with predictions of continuum theory. Simulations 
are at a a = 0.045 and various L y — L z . 
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